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Abstract 

Complex network dynamics have been analyzed with models of systems of coupled switches or systems of 
coupled oscillators. However, many complex systems are composed of components with diverse dynamics 
whose interactions drive the system's evolution. We, therefore, introduce a new modeling framework that 
describes the dynamics of networks composed of both oscillators and switches. Both oscillator synchro- 
nization and switch stability are preserved in these heterogeneous, coupled networks. Furthermore, this 
model recapitulates the qualitative dynamics for the yeast cell cycle consistent with the hypothesized 
dynamics resulting from decomposition of the regulatory network into dynamic motifs. Introducing feed- 
back into the cell-cycle network induces qualitative dynamics analogous to limitless replicative potential 
that is a hallmark of cancer. As a result, the proposed model of switch and oscillator coupling provides 
the ability to incorporate mechanisms that underlie the synchronized stimulus response ubiquitous in 
biochemical systems. 

Introduction 

The dynamics in systems ranging from intercellular gene regulation to organogenesis are driven by com- 
plex interactions (represented as edges) in subcomponents (represented as nodes) in networks. If the 
structure of these networks is known, network-wide models of coupled systems have been applied to pre- 
dict their qualitative dynamics. For example, models of coupled switches based upon Glass networks [l] 
have been applied to model systems such as neuronal synapses \2 and gene regulatory networks fsl. 
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Similarly, models of coupled oscillators along networks based upon the Kuramoto model [4] have been 
used to model synchronization of oscillators in diverse systems reviewed in |5l. In biochemical systems, in 
vivo oscillator synchronization has been observed in synthetic oscillatory fluorescent bacteria (6][7|, yeast 
gene regulatory networks [8||9], and human cell fate decisions llOj. Such spontaneous synchronization has 



also been attributed to the development of the mammalian cardiac pacemaker cells (reviewed in 11 ) and 
cortical systems (reviewed in !ll2|) including notably the circadian pacemaker (e.g., [13| ). More recently, 
these network models have been found to be insufficient to model more complex dynamics in neuronal 



information transfer 12 14 17 and cardiac arrhythmias [18-21 . These limitations extend to physical 



systems, such as the coupled lasers studied in 22 . Therefore, numerous studies have modified these 



network models to account for evolving networks 15 23-28 , dynamic frequencies 15 29 30 , or phase 



delays 16 31-33 . However, these mathematical modifications typically do not encode the mechanism 
underlying the limitations in the Kuramoto and Glass network models. 

We hypothesize that the observed limitations in the standard Kuramoto and Glass models arise 
from their exclusion of coupling components with qualitatively different dynamics. Several studies have 
inferred that biochemical systems contain "network motifs" with both oscillatory and switch-like dynamics 



34 35 . The dynamics of these motifs are inferred from the topology of subgraphs in the networks of 



these systems. Their structures are statistically overrepresented in biochemical networks [36j[37] such as 
intracellular regulatory networks [38 , implicating evolutionary preservation (and thus utility) of these 



network motifs | 39j. The dynamics of these motifs have been used to model yeast cell cycle regulation 40 



and have been further confirmed in synthetic, designed biochemical circuits (reviewed in [41]). Because 
these heterogeneous network motifs are all identified as components within a single biochemical network. 



their interactions must drive the global dynamics of the network 42 . Previously, [43[ have shown that 
coupling small sets of heterogeneous network motifs ensures the robustness of motif dynamics and |42j have 
shown that coupling networks changes their dynamics in isolation. However, the network-level dynamics 
that result from coupling oscillatory and switch-like components have not been studied comprehensively. 

In this paper, we develop a theoretical framework to quantify the network-wide dynamics resulting 
from coupling switches and oscillators. This model is based upon introducing cross-coupling between 
the Kuramoto and Glass models, due to their wide success in modeling the dynamics in networks of 
oscillators and networks of switches, respectively. Simulations with the proposed model across state- 
space in an all-to-all network yields four operational states: (1) switches remain "on" and oscillators 
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synchronize, (2) switches are "off" and osciUators freeze, (3) switches fluctuate in sync with oscillators, 
and (4) switches fluctuate transitionally until oscillators freeze. Further application of our model to the 
network motifs identified in yeast in \4A^ recapitulates the qualitative dynamics of the system observed 
in that study. However, a simple rewiring of this cell-cycle network that introduces feedback causes a 
cancer-like sustained re-activation of the cell cycle machinery without regard for external signal growth 
signals. These dynamics suggest that modeling cross-motif coupling may predict critical processes in the 
dynamics of biochemical networks with minimal parameterization. 

The Kuramoto model of coupled oscillators 

Quantitative studies of coupled oscillators often apply the Kuramoto model of AI oscillators coupled in 
an all-to-all network. In this model, the change in time 9i of the phase of the i*^ oscillator, is governed 

by 

M 

9.^6j.+ '^J2'in{9,~0,), (1) 

where loi is the natural frequency of the oscillator and Kg^g > is the coupling strength of the 
oscillators [4j. Typically, the uji values are drawn from a normal distribution centered at with variance 

In the Kuramoto model, the phases of the oscillators will synchronize if Kg g is above a threshold 
coupling strength Kg g. Such synchronization is quantified with the mean field of the oscillators as 

M 

Here ij) is the average phase of the oscillators and the coherence rg represents the spread of the oscillators 
from that average phase. Based upon eq. ([2]), = 1 if each 9i — ip and = if the values of 9i are 
distributed miiformly between [0, 27r) [45] . 
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Glass networks of coupled switches 

Coupled sets of N switches, which adopt one of a set of binary states, are modeled with Glass networks [l] . 
These models describe the evolution of the switch (xi) as follows 

Xi = -Xi + F,{xi,X2, ■ ■ ■ ,Xn) , a.Tid (3) 

ii = if Xi < 0; 1 otherwise, (4) 

where Xi represents the change in time of the value of each Xi , which are unobservable continuous variables 
that control the time of switching between observable, discrete states in Xi. In this model, Fi describes 
the change in state of the i^^ switch due to the coupling with the other N switches in the network [ij. 
In specified network structures and functions Fi, such Glass networks can exhibit complex dynamics, 
including periodic and aperiodic orbits (e.g., [46| ) . 

One type of Glass network, called a Hopfield network , has dynamics applicable to the smooth-decay 
of signal in biochemical switches The Hopfield model lets 

AT 
J = l 

where Wij takes values between —1 and 1 representing the relative strength of the connection between 
switches i and j, k^^x is the magnitude of coupling strengths, and t,; the threshold for switch activation. 
Similar to the Kuramoto model, sets of the switches will synchronize for Kx^x above a threshold kx.x in 
appropriate network topologies. 

Results 

Network model of coupled oscillators and switches 

By combining the established models for switches and oscillators, we model the dynamics of the hetero- 
geneous system of coupled switches and oscillators in systems including biochemical networks with the 
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following set of equations: 

Xi = -Xi+Gi{xi,X2,---,XN,0i,92,...,9M) (6) 

01 = oji{xi,X2,...,xn) (7) 

+ Hi {xi,X2, . ■ .,XM,9i,02, ■ ■ ..Om) ■ 

Here, eq. (|6| is analogous to the Glass network in eq. ([s]) and Xi is defined according to eq. Q. 

In this study, we explore a case of the switch-oscillator model in cqs. Q and ([t]) which contains an 
all-to-all network that couples the Kuramoto model, eq. ([l]), and Hopfield network, eqs. (|3| - ([s]), as 
follows 

N M 

= -h ^ ^ ^ ^ e'fc - n, (8) 



N ^ ' M 

j^. I k=l 
M 



K0, 



M 

k=l 



J2smiek'9i), (9) 



N 

^(ijw, - w,), (10) 



where n^ g and ng^^ si's cross-component coupling strengths. In eq. (10), is the time- varying frequency 



of the l^^^ oscillator resulting from switch coupling, with initial values uji {t — 0) — uii, and 

1 if < 61/ < TT 

n = s ■ (11) 

otherwise 

In this system, zero values of the cross-coupling parameters Kx,e and Kg^x cause the model to reduce to 
the standard uncoupled Kuramoto and Hopfield models. Similar decoupling of the models occurs if the 
switch and oscillator systems are at vastly different timescales, determined by the and uJi parameters. 



respectively. The transformation in eq. (Ill facilitates comparable switch-like dynamics in the oscillators 
when they interact with switches in eq. ([8|. Nonzero switch-oscillator (nx.e) interactions will cause an 
oscillator in the "up" part {9i — 1) of its cycle to feed energy into the switch in question, nudging it 
towards the "on" state if off or delaying its decay if already on. Similarly, an "on" switch with a nonzero 
oscillator-switch interaction (Kg^^) will feed energy into the oscillators causing them to cycle at their 
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natural frequency if coupled to that switch. 

By thus incorporating coupling between switches and oscillators within the framework established by 



the standard Kuramoto and Hopfield models, the dynamics of our model in cqs. ([8|) - (10) can be analyzed 
within the framework of these well established models. Similar to analysis of the Kuramoto model and 
Glass network, we summarize the dynamics of our system using order parameters. For oscillators, we 
utilize the order parameter defined in eq. ([2|. We introduce a new order parameter 

^"W-]?7E^ (12) 

Ad ^-^ UJu 

that tracks how closely each individual oscillator's frequency lj^ corresponds to the natural frequency uj^. 
Analogously, we measure the fraction of switches that are in the "on" position at a given time using a 
switch-switch order parameter defined by 



Both of these functions will have a maximum of 1 when all switches are on, and minimum if all switches 
are off. 

Simulation results in all-to-all networks 

We first explore the qualitative dynamics of the heterogeneous system through numerical simulations 
in all-to-all networks. We limited these simulations to all-to-all networks, because of the ability of this 
network topology to describe the qualitative dynamics from the Kuramoto model. These simulations 
explore the majority of parameter space defined by Kx,xi i^x,e, f^e.e, Kg,x, and cr^. Specifically, we select 
Kx^x = 1 < Kx,x to ensure that switches are able to turn off without appropriate stimulation from the 
oscillators. We consider the effects of switches on oscillators for values of Kg,g both above and below 
the Kuramoto threshold kg^g. Figure [T] plots the time-dependent order parameters observed in the four 



qualitative states observed in simulations of the coupled model eqs. ([8|) - ( 10 ) that are reflective of the 
qualitative dynamics observed in simulations with these parameter values. Supplemental videos S1-S4 
further summarize the results of these simulations. We note that these four states were the only qualitative 
states observed for our coupled model in all-to-all networks simulated according to the description in the 
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Methods section. Because r, n^fi^ and a^^ all control the relative timing of switches and oscillators, their 
values were selected in these simulations to optimize visualization in the supplemental videos. When 
exploring the effect of timing on the system dynamics, we hold r and n^fi fixed while varying a^. 
Figure [2] shows the probability of observing the states in Figures l(a)||f (c) in 100 simulations of all-to-all 



systems containing 100 switches and oscillators as a function of kq^x and Ui^. Because of their common 
control of system timing, we would obtain comparable distributions when varying either r or n^fi instead 

of (T„. 

The coupled system preserves synchronization in both oscillators and switches. 



Figure 1(a) shows a state of the model in which the switches are all in the "on" state and oscillators 
are synchronized {r^ near 1, rg near 1, and ^ oscillating between [0, 27r) periodically). While such 
synchronization is observed in the uncoupled Hopfield and Kuramoto models, the oscillator-switch cross 
coupling extends the region of parameter space over which this synchronization occurs. Specifically, 
modest values of n^.e can induce sustained switch activity for parameter values of Kx.x in which switches 
would decay in the uncoupled system. Furthermore, this switch synchronization will occur for all values 
of Kx,x in which synchronization occurs in the uncoupled Hopfield model (i.e., all Kx^x larger than a 
threshold value kx^x) because the oscillators only contribute positively to the derivative in eq. Q in our 
model. On the other hand, no value of kq^x will cause oscillator phases to synchronize if Kg^ is below the 
critical coupling parameter for the pure Kuramoto model [kg_g). However, there are parameter regimes 
in which this synchronization occurs stochastically, depending on the initial values selected for Xi, 0j, 



and (jjj (Figure 2(a)). In these cases, the average decrease in oscillator natural frequencies caused by 
decreasing Kg x or cr^ will increase the effective period of oscillators, thereby increasing the probability of 
switches being locked in the "on" state and oscillator synchronization in the heterogeneous system. 

Coupling switches to unsynchronized oscillators can freeze network-wide dynamics. 



Figure 1(b) depicts a model state in which switches are all "off' (r^, near zero) and oscillators "freeze": 
each 9j (t) = (t) = ^ for some constant values for all t beyond the preliminary freezing time t y . While 
the decaying switches are observed in an uncoupled Hopfield model, the freezing oscillators cannot be 
simulated in the uncoupled Kuramoto model. Such oscillator freezing will occur whenever the oscillators 



decay to the "off' state by virtue of the coupling of the oscillators to switches through the ujj in eq. ( 10 ) 
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Specifically, this frozen state can occur whenever k^^^ < Kx^x depending on the values of Xi, 0j, and 
ujj. However, the probability of selecting these initial states is decreased when the heterogeneity of the 
oscillators increases through incomplete synchronization {rg{t) < 1) or increased a^o (Figure [2(b)[ ). In 
these cases, a single oscillator in the "up" phase {6 = 1) can contribute positively to the switch states, 
forcing the system out of this frozen state. The probability of obtaining this frozen model state further 
depends on the relative timing of switch decay and oscillator freezing. Specifically, the probability of 
obtaining the frozen state decreases with the average oscillator frequency, determined predominantly by 



the parameter Kg,x (Figure 2(b) ) 



Coupling switches to synchronized oscillators can induce synchronized oscillations in switches. 

An additional consequence of coupling switches and oscillators in a state in which switches vacillate 



between all "on" and all "off' along with the synchronized oscillator frequency (Figure 1(c) I. This 
oscillatory synchronization occurs when the pure Hopfield model would turn switches "ofF {kx^x < 'ix,x), 
the pure Kuramoto model would induce oscillator synchronization {ng^g > kg^g), and the timing between 
the oscillators and switches are balanced such that the average period of the coupled oscillators is slightly 



less than the average decay time of the system of switches. Figure 2(c) shows that this balance in switch- 
oscillator timescales increases with decreasing ng^x and depends non-monotonically on ct^j. As we see in 
the plot of (t) in Figure |l(c) the average oscillator natural frequencies will decrease towards the end 



of the "down" phase in response to switches turning off, and then increase to their full natural values in 
the "up" phase as switches turn back "on". Therefore, if synchronized oscillator period is too slow (i.e.. 



(7i^ is too large), the system will tend to be locked in the "on" state (Figure 2(a) ); if too fast (i.e., too 



small) the system will tend to be locked in the "off' state (Figure 2(b) ) 



Synchronization of network-wide oscillations may be transitory. 

Oscillatory behavior in the switches is also observed for unsynchronized oscillators {Kg^g < kg^g) as 



depicted in Figure 1(d) In this case, the value of Kg^x must be large enough to enable switches to freeze 
the oscillators' phases. However, because the oscillators are uncoupled, a small subset of oscillators in the 
"up" phase can drive the switches to turn on for large-enough values of Kx^g. These switch oscillations are 
transitory, ending when at last the switch coupling dominates the system and induces all of the oscillators 



to freeze. For unsynchronized oscillators in the parameter range of Figure 1(d) the transitional oscillations 
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in the switch state occurs regularly in 21 of 100 simulations. In 8 of these 21 simulations, the switch 
state turns "on" after decaying at least twice. More rarely, transitory changes in switch state may be 
induced by a similar mechanism in simulations for which Kg^g > ke^e and switches ultimately settle on 
the all "on" or all "ofF states. 

System size affects the distribution of qualitative dynamics 

We also explored the dynamics of the coupled system for networks of sizes ranging from N = M = 10 to 
N — M = 500 nodes, described in the methods. For networks of all sizes, we observe that the dynamics 
of the system was limited to the four qualitative behaviors observed for networks of size N — M ^ 100 
depicted in Figure [ij However, the system size does have a notable effect on the frequency with which 
each of these behaviors occurs. Supplemental Figures S5-S7 plot the observed frequencies for each of the 
network sizes as a function of the ng,x and cr^ values considered in Figure [2] 

When Kg^x = 0.01, the observed frequencies of the system states depend most strongly on network 
size in simulations using the smallest value oi = \ is also small (Supplemental Figure S5). In this 
case, the probability of observing the system with synchronized oscillatory dynamics in both switches and 
oscillators decays as the network grows. Both the state in which the switches are on and oscillators are 
synchronized and the state in which the switches are off and oscillators are frozen have with compensatory 
increases in probability (Figure |3]). The relative probability of obtaining the frozen state increases, with 
notable decay in the probability of obtaining the state in which switches are "on" and oscillators are 
synchronized in large networks. 

On the other hand, when k^.x — 1, the system size has the greatest influence on the resulting dynamics 
for large values of a^^ (Supplemental Figure S7). In this case, the system changes from containing mostly 
switches in the on state and synchronized oscillators to switches that are entirely in the "frozen" state for 
large network sizes (Figure [4|. We hypothesize that the system is forced into the frozen state in larger 
networks because of increased oscillator synchronization in large networks. Therefore, small networks 
would have a higher probability of having few oscillators that are unsynchronized and in the "up" phase 
{9 = 1), causing the switches to turn "on" [x = 1) due to the structure of eq. ([s]) as was discussed 
previously. Furthermore, the rare oscillations observed in both switches and oscillators when ng^^ — 1 
occur only when the network is small. Intermediate values of Kg^x = 0.1 show similar changes to those 
described for Kg^x = 0.01 when a^^j — 1 and to those described for Kg^x = 1 when CTi^ = 10 (Supplemental 
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Figure S6). 



The heterogeneous network models quahtative dynamics of the yeast cell cycle 
derived from network motifs. 



Previous work by 47 make the ceh cycle processes controlling mitotic division of fission yeast Schizosac- 
charomyces pombe cells provides an optimal system in which to apply our model. The biochemical 
reactions responsible for driving the cell cycle are well understood and the resulting dynamics in each 



of the stages of the cell cycle have been characterized extensively in 40 44 47 . The cell cycle machin- 
ery in mitosis is divided into four, sequential stages: phase 1 is a gap or rest phase (Gl); phase 2 is a 
DNA synthesis stage (S); phase 3 is an additional gap stage (G2); and phase 4 is the mitotic division 
stage (M). Previously, [47| observed that the dynamics of the yeast cell cycle can be divided into three 
sequentially interacting modules, triggered by a signal based upon cell size: (1) Gl/S transitions with 
a toggle-switch, (2) S/G2 transitions with a toggle-switch, and (3) G2/M transitions with an oscillator. 
Although the specific timing differs from f47], we observe similar qualitative dynamics to that observed 
in 47 when applying our heterogeneous model to evolve the state of these cell cycle stages (Figure [5]) 
as described in the Methods section. We note that the response in this system is consistent with the 



transitory oscillations observed in Figure 1(d) in the case of all-to-all coupling. We also modeled this 
cell-cycle system in a rewired-network, in which the G2/M transitions feedback into Gl/S (Figure [6]). 
In this case, we observe sustained reactivation of the cell cycle regardless of the external signal. These 



dynamics are analogous to the synchronized dynamics in Figure 2(c) and consistent with cell growth 
arising from re- wiring biochemical reactions in cancer cells [48j. 



Discussion 

Our model of coupled switches and oscillators in all-to-all networks demonstrates that networks with 
components having heterogeneous dynamics can exhibit synchronization similar to that observed in ho- 
mogeneous systems. As is the case in homogeneous models (e.g., [49}j52] ), we expect analogous synchro- 
nization to hold in small- world, biochemical network topologies (e.g., [53| ). However, these alternative 
topologies would likely change the probability of observing each of the qualitative model behaviors similar 
to the observed dependence of probabilities in network size. In this alternative network topologies, the 
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qualitative states of the network model may have greater variability in small network sizes in accordance 
with the findings of 54|. Finally, in these topologies the heterogeneous model could yield additional, 
complex qualitative dynamic states, resulting from the complex dynamics that they cause in models of 



coupled switches alone 46 



While uncoupled network motifs may adopt switch-like or oscillatory dynamics, coupling between 
these components can induce switch-like behavior in oscillators and oscillatory behavior in switches. 
These qualitative changes in component dynamics occur stochastically, depending on the distribution of 
frequencies and switch states. They are more likely to occur in simulations with an imbalance in relative 
timescales, in which the dynamics of the faster network motif will dominate the system. Similarly, when 
Ke,x and (7^^ are both small, the coordinated oscillations in the switches and oscillators that occur in 
frequently small networks are largely eliminated in larger networks. We hypothesize that this larger 
network effectively increases the range of natural frequencies and phases, making the simulation less 
likely to have the constrained distribution required to obtain such synchronized oscillations. We can 
expect that biological systems have evolved components according to these distributions to ensure the 
robustness of the dynamics in the system. For example, multiple proteins can often serve similar functions 
in cell signaling pathways, which would increase the system size and decrease the probability of transient 
behaviors in our model. This robustness will be further ensured through the sheer size of most biochemical 
systems. For example, in humans yeast two-hybrid maps and metabolic network maps both contain on 



the order of thousands of interactions between thousands of species 53 



Furthermore, we have also observed that the heterogeneous network model will freeze the oscillator 
dynamics in the presence of inactive switches and then subsequently activate in synchrony in the presence 
of active switches. As a result, our model provides a natural mechanism for the coordination of complex 
machinery such as the initiation of cell-cycle dynamics. For example, when we apply our model to the 



yeast cell cycle motifs in 47 , we recapitulate the qualitative dynamics of delayed initiation of stages 



of the cell cycle observed in simulations with differential equations of the regulatory dynamics in 47 



Additional tuning of the model parameters or incorporation of additional cell cycle checkpoints would 



facilitate a precise match of the timing of 47 . Because parameters are defined for modules and their 
interaction, our model requires far fewer rate parameters than any differential equation model of sets 
of biochemical reactions of the yeast cell cycle. Generally, the oscillator in the final G2/M step of the 
cell cycle is active only when the series of switches in the previous steps of the cell cycle are activated. 
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consistent with the transient dynamics observed in our network model. However, rewiring the network 
to introduce feedback from the G2/M stage to the Gl/S stage of the ceh cycle will cause the modeled 
cell cycle machinery to engage continually without regard to the external growth signals, consistent with 



the malignant rewring in cancer cells 48 . Similar to the oscillatory behavior induced in switches in 
simulations in all-to-all networks, this small modification to the topology of cell cycle interactions altered 
the resulting dynamics of the network motifs for the Gl/S and S/G2 motifs. We, therefore, hypothesize 
that motif dynamics predicted by the structure of subgraphs may not accurately describe their in vivo 



dynamics if considered in isolation, consistent with the hypothesis in 55 and findings of 42 



We observed that the switches in the cell cycle block activation of the yeast cell cycle when no external 
signal is present. Similarly, when part of the larger but sparse networks that compose biochemical systems 



53 , inactive switches would effectively destroy links between nodes on the network. As a result, the 



proposed heterogeneous model provides a potential mechanism for Kuramoto-based models with evolving 



network topologies such as [15 23 -28 . Similarly, we observed that the intermediate switches delay the 
oscillations in the final G2/M motif in the simulated yeast cell cycle. As a result, we hypothesize that 
coupling switches to oscillators through their frequencies in this model also provides a natural mechanism 



for extensions of the Kuramoto model with dynamic frequencies 15 29 30 or phase delays 16 31 - 33 



The heterogeneous network model described in this paper facilitates characterization of the dynam- 
ics of complex, biochemical systems by abstracting the dynamics of their composite motifs such as the 
yeast cell cycle based upon [4^. We note that the proposed heterogeneous network model is determin- 
istic once the initial values of all the switches and oscillator frequencies have been specified. However, 
many intracellular reactions (e.g., [56| ) and neuronal systems (reviewed in [57[|58[ ) evolve stochastically. 
In these cases, the Hopfield networks used to model the switches could be replaced with probabilistic 
Boolean networks [3] and the oscillators evolved with stochastic solvers such as the stochastic simulation 



algorithm (reviewed in 59 ), integrated with the methodology developed in |60|. Similar modifications 
could also extend the heterogeneous model to incorporate coupling with components of additional dy- 



namics pertinent to biochemical systems, such as those of the network motifs enumerated in f34 35 44 



These studies would also ideally consider the dynamics of the heterogeneous network model in additional 
small- world and random network topologies, as well as the topologies defined by neuronal systems and 
gene regulatory networks. 
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Materials and Methods 

Numerical simulations in the all-to-all network 

In this study, we analyze the range of possible dynamics of the coupled, heterogeneous networks by 
applying this model to all-to-all networks. Analyses were performed for networks with equal number 
of switches and oscillators {N — M) of sizes 10, 50, 100, 200, and 500. All simulations are run one 
hundred times from random initial conditions for the state of switches {xi^i = 1, . . . ,N) and oscillators 
{Oj,j = 1, . . . , Af), drawn from a Gaussian distribution and a uniform distribution on [0, 27r), respectively. 
Similarly, oscillator natural frequencies are drawn randomly from a Gaussian distribution of mean zero 
and standard deviation parameter a^^. Simulations of 100 seconds (in the arbitrary units of the model), 
with a time step of 0.01 seconds were found sufficient to reflect the range of possible model behaviors and 
verify consistency across initial conditions. The behavior of each simulation is summarized based on the 
time-dependent order parameters re (t) and ip (t), (t), and Tx (t). 



Numerical simulations of the yeast cell cycle 



Based upon 47 , we model the yeast cell cycle as an initiating external signal (namely the cell size), 
coupled to a toggle switch representing the transition between Gl/S, a toggle switch representing the 
transition between S/G2, and an oscillator representing the transition from G2/M. While the external 
signal is incorporated into the model with coupling to the other switches in eq. (|6| , its state is not updated 
by the model. The duration of this external signal is set at 10 simulated minutes, based upon [47^ . 
Similarly, the initial values of the hidden variable x for the switches in the Gl/S and S/G2 modules are 
set at -0.5, t to 1, and k^^x to 2 to reproduce the approximate 10 minute duration of these switches in [47). 
The natural frequency is for the G2/M module set to |^min^^ to likewise reflect the timescale reported 
while the remainder of the coupling parameters are left untuned, set to Ke,x = f^x.e = i^e.o = 2 



47 



because we sought only to reproduce the qualitative dynamics of the 47 model. The rewiring in the 
system with enduring cell cycle activation is achieved by adding an edge from the module for G2/M to 
the switch in the Gl/S module. 
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Summary of the qualitative dynamics of the heterogeneous network model of 
eqs. ^ - ([To]). In all fi gures, top-panel shows temporal evolution of the mean field 



statistics {rg black, solid; green, dashed; and blue, dash-dotted) and the bottom- 
panel shows the evolution of the mean phase tp (red, solid), (a) Oscillators synchronize 
and switches stay "on" {k^^x — H, n^.e — 1.5, Kg^x — 1, — 40, and a^j ~ 10), 
(b) oscillators freeze (as evidenced by unchanging ip) and switches stay "off' {kx,x = 1, 
Kx,e — 1-5, Ke^x = Ij ng^e = 40, and — 10), (c) oscillators synchronize and switches 
oscillate {kx^x = 1, Kx.e = 160, Kg^x = 0.2, Kg^g = 42, and auj = 3), and (d) transitory 
oscillations in oscillators and switches (kx^x =0.1, Kx.e — 1.4, Kg_x — 2, Kg^g = 1.8, and 

<Ju = lO) 

Percentage of simulations in which the qualitative dynamics in Figure [T] occur. 

In (a) oscillators synchronize and switches are "on" , in (b) oscillators freeze and switches 
are "off', and in (c) switches vary with oscillators vs a^^ for Kg^x = 0.01 (solid), Kg,x = 0.1 

(dotted) and Kg,x — 1 (dashed). Kx^x = 1 < f^x^x, i^x,0 — 1-5, and Kg_g = 40 > kg^ 

Dependence on network size for qualitative states for Kg,x — 0.01 and a^^ = 
1. Percentage of simulations in which the qualitative dynamics have switches off and 
oscillators frozen (blue, solid), switches on and oscillators synchronized (green, dashed), 

and oscillatory switches and synchronized oscillators (red, dotted) 

Dependence on network size for qualitative states for ng^x = 1 and = 10. 
Percentage of simulations with qualitative dynamics plotted as described in Figure |3] . . . 
Simulated dynamics of the yeast cell cycle Evolution of the states of the cell cycle 
modules (Gl/S top, green dashed; S/G2 top, red dotted; G2/M bottom, black) in response 

to an external stimulus to initiate the cell cycle (top, blue solid) 

Simulated dynamics of an aberrant cell cycle network. As for Figure [5] with a 
network topology linking the G2/M module to the Gl/S module 
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Figure 1. Summary of the qualitative dynamics of the heterogeneous network model of 
eqs. @ - In all fi gures, top-panel shows temporal evolution of the mean field statistics (rg 

black, solid; green, dashed; and r^^ blue, dash-dotted) and the bottom-panel shows the evolution of 
the mean phase ip (red, solid), (a) Oscillators synchronize and switches stay "on" {kx,x = H, fix,e = 1-5, 
K-d.x = 1, ng,e = 40, and cr^^ = 10), (b) oscillators freeze (as evidenced by unchanging ^) and switches 
stay "ofF {kx,x = 1, Hx,9 = 1-5, ng^^ = 1, = 40, and cr^j = 10), (c) oscillators synchronize and 
switches oscillate {kx,x = 1, ^^^e = 160, Kg = 0.2, Kg e = 42, and a^^ =3), and (d) transitory 
oscillations in oscillators and switches {kx,x = 0-1, = 1.4, Kg ^ — 2, Kg g — 1.8, and — 10). 
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(a) (b) (c) 



Figure 2. Percentage of simulations in which the qualitative dynamics in Figure [ij occur. 

In (a) oscillators synchronize and switches are "on" , in (b) oscillators freeze and switches are "off' , and 
in (c) switches vary with oscillators vs a^^ for kq^x ~ 0.01 (solid), ng^x ~ 0.1 (dotted) and ng^x — 1 
(dashed). Kx,x = 1 < kx,x, Kx,e = 1-5, and Ke,e = 40 > kg^g. 

Supplemental Figure Captions 

Figure S5. Dependence of dynamics on network size for Hig,x — 0.01. Number of simulations 
(of 100) for which the switches are off and oscillators are frozen (left panel), the switches are on and 
the oscillators are synchronized (center panel), and both the oscillators and switches have synchronized 
oscillations (right). 

Figure S6. Dependence of dynamics on network size for kq x = 0.1. As for Supplemental 
Figure S5. 

Figure S7. Dependence of dynamics on network size for ng x — 1- As for Supplemental 
Figure S5. 
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Figure 3. Dependence on network size for qualitative states for Kg, a; = 0.01 and ct^ = 1. 

Percentage of simulations in which the qualitative dynamics have switches off and oscillators frozen 
(blue, solid), switches on and oscillators synchronized (green, dashed), and oscillatory switches and 
synchronized oscillators (red, dotted). 
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Figure 4. Dependence on network size for qualitative states for Kg^x = 1 and = 10. 
Percentage of simulations with qualitative dynamics plotted as described in Figure [Sl 
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Figure 5. Simulated dynamics of the yeast cell cycle Evolution of the states of the ceU cycle 
modules (Gl/S top, green dashed; S/G2 top, red dotted; G2/M bottom, black) in response to an 
external stimulus to initiate the cell cycle (top, blue solid) 




10 20 30 40 50 60 70 80 90 100 

time (min) 



Figure 6. Simulated dynamics of an aberrant cell cycle network. As for Figure [5| with a 
network topology linking the G2/M module to the Gl/S module. 



